Ejemplo 7

EDA

milk <- read.table("milk.txt",header=TRUE)
milk$t<-1970:1982
n<-nrow(milk)
or<-order(milk$x)
plot(milk$x[or],milk$y[or],type="l")
text(milk$x[or],milk$y[or],labels=milk$t[or],cex=0.5,col=2)

plot(milk$t,milk$y,type="l")

plot(milk$t,milk$x,type="l")

Modelo O

#-Defining data-
m<-2
# data<-list("n"=n,"m"=m,"y"=milk$y,"x"=milk$x,"t"=milk$t)
# data<-list("n"=n,"m"=m,"y"=milk$y,"x"=milk$x/max(milk$x),"t"=milk$t/max(milk$t))
data<-list("n"=n,"m"=m,"y"=scale(milk$y)[1:n],"x"=scale(milk$x)[1:n],"t"=scale(milk$t)[1:n])
# data<-list("n"=n,"m"=m,"y"=c(scale(milk$y)[1:(n-2)],NA,NA),"x"=scale(milk$x)[1:n],"t"=scale(milk$t)[1:n])

#-Defining inits-
inits<-function(){list(alpha=0,beta=0,tau=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,5),tau=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,n+m),tau.y=1,tau.b=1,yf1=rep(0,n+m))}
# inits<-function(){list(beta=rep(0,n+m),tau.y=1,yf1=rep(0,n+m))}
# inits<-function(){list(alpha=rep(0,n+m),beta=rep(0,n+m),tau.y=1,tau.b=1,tau.a=1,yf1=rep(0,n+m))}
# inits<-function(){list(beta=rep(0,n+m),tau.y=1,tau.b=1,yf1=rep(0,n+m),g=0)}
# inits<-function(){list(beta=rep(0,n),tau.y=1,tau.b=1,yf1=rep(0,n),g=1)}
# inits<-function(){list(beta=rep(0,n),tau.y=1,tau.b=1,yf1=rep(0,n))}
# inits<-function(){list(beta=rep(0,n),yf1=rep(0,n))}

#-Selecting parameters to monitor-
parameters<-c("beta","tau","yf1")
# parameters<-c("beta","tau.y","tau.b","yf1","g")
# parameters<-c("beta","tau.y","tau.b","yf1")
# parameters<-c("alpha","beta","tau.y","tau.b","tau.a","yf1")
# #parameters<-c("beta","tau.y","tau.b","yf1")
# parameters<-c("beta","yf1")

#-Running code-
#OpenBUGS
ej7o.bugs<-bugs(data,inits,parameters,model.file="Ej7o.txt",
               n.iter=50000,n.chains=2,n.burnin=5000)
# ej7a.bugs<-bugs(data,inits,parameters,model.file="Ej7a.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000)
# ej7b.bugs<-bugs(data,inits,parameters,model.file="Ej7b.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000)
# ej7c.bugs<-bugs(data,inits,parameters,model.file="Ej7c.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000,debug=TRUE)
# ej7d.bugs<-bugs(data,inits,parameters,model.file="Ej7d.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000)
#JAGS
ej7o.jags<-jags(data,inits,parameters,model.file="Ej7o.txt",
               n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)

## module glm loaded

## Warning in jags.model(model.file, data = data, inits = init.values,
## n.chains = n.chains, : Unused variable "m" in data

## Warning in jags.model(model.file, data = data, inits = init.values,
## n.chains = n.chains, : Unused variable "t" in data

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 13
##    Unobserved stochastic nodes: 16
##    Total graph size: 69
## 
## Initializing model

# ej7a.jags<-jags(data,inits,parameters,model.file="Ej7a.txt",
#               n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)
# ej7b.jags<-jags(data,inits,parameters,model.file="Ej7b.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)
# ej7c.jags<-jags(data,inits,parameters,model.file="Ej7c.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)

exploracionMCMC <- function(modelo){
  if("BUGSoutput" %in% names(modelo)){
    modelo <- modelo$BUGSoutput
  }
  out <- modelo$sims.list
  z <- out$beta[]
  par(mfrow=c(2,2))
  plot(z,type="l")
  plot(cumsum(z)/(1:length(z)),type="l")
  hist(z,freq=FALSE)
  acf(z)
  par(mfrow=c(2,2))
}
exploracionMCMC(ej7o.bugs)

#-Monitoring chain-

#Traza de la cadena
# traceplot(ej7o.bugs)

dicMCMC(ej7o.bugs)

## [1] 38.68

dicMCMC(ej7o.jags)

## [1] 39.26374

ej7o.bugs$summary[1:2,]

##            mean        sd    2.5%     25%      50%       75%      97.5%
## beta -0.5045959 0.2886381 -1.0810 -0.6862 -0.50445 -0.322575 0.06927025
## tau   1.2281515 0.5257877  0.4243  0.8455  1.15600  1.531000 2.44502500
##          Rhat n.eff
## beta 1.000987 90000
## tau  1.000988 90000

#Resumen (estimadores)
#OpenBUGS
out.sum<-ej7o.bugs$summary

#JAGS
# out.sum<-ej7a.sim$BUGSoutput$summary

#Tabla resumen
out.sum.t<-out.sum[grep("beta",rownames(out.sum)),c(1,3,7)]
# out.sum.t<-cbind(out.sum.t,apply(out$beta,2,prob))
# dimnames(out.sum.t)[[2]][4]<-"prob"
print(out.sum.t)

##        mean        2.5%       97.5% 
## -0.50459588 -1.08100000  0.06927025

#DIC
out.dic<-ej7o.bugs$DIC
# out.dic<-ej7a.sim$BUGSoutput$DIC
print(out.dic)

## [1] 38.68

#Predictions
out.yf<-out.sum[grep("yf1",rownames(out.sum)),]
ymin<-min(data$y,out.yf[,c(1,3,7)])
ymax<-max(data$y,out.yf[,c(1,3,7)])
xmin<-min(data$t)
xmax<-max(data$t+m)

#x vs. y
par(mfrow=c(1,1))
plot(data$x,data$y,type="p",col="grey50",ylim=c(ymin,ymax))
points(data$x,out.yf[,1],col=2,pch=16,cex=0.5)
segments(data$x,out.yf[,3],data$x,out.yf[,7],col=2)

#t vs y
par(mfrow=c(1,1))
plot(data$t,data$y,type="b",col="grey80",ylim=c(ymin,ymax),xlim=c(xmin,xmax))
lines(data$t,out.yf[1:n,1],col=2)
lines(data$t,out.yf[1:n,3],col=2,lty=2)
lines(data$t,out.yf[1:n,7],col=2,lty=2)

# lines(data$t[n]:(data$t[n]+m),out.yf[n:(n+m),1],col=4)
# lines(data$t[n]:(data$t[n]+m),out.yf[n:(n+m),3],col=4,lty=2)
# lines(data$t[n]:(data$t[n]+m),out.yf[n:(n+m),7],col=4,lty=2)
#betas
# out.beta<-out.sum[grep("beta",rownames(out.sum)),]
# ymin<-min(out.beta[,c(1,3,7)])
# ymax<-max(out.beta[,c(1,3,7)])
# plot(out.beta[,1],type="l",ylim=c(ymin,ymax))
# lines(out.beta[,3],lty=2)
# lines(out.beta[,7],lty=2)

Modelo A

#-Defining data-
m<-2
# data<-list("n"=n,"m"=m,"y"=milk$y,"x"=milk$x,"t"=milk$t)
# data<-list("n"=n,"m"=m,"y"=milk$y,"x"=milk$x/max(milk$x),"t"=milk$t/max(milk$t))
data<-list("n"=n,"m"=m,"y"=scale(milk$y)[1:n],"x"=scale(milk$x)[1:n],"t"=scale(milk$t)[1:n])
# data<-list("n"=n,"m"=m,"y"=c(scale(milk$y)[1:(n-2)],NA,NA),"x"=scale(milk$x)[1:n],"t"=scale(milk$t)[1:n])

#-Defining inits-
# inits<-function(){list(alpha=0,beta=0,tau=1,yf1=rep(1,n))}
inits<-function(){list(beta=rep(0,5),tau=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,n+m),tau.y=1,tau.b=1,yf1=rep(0,n+m))}
# inits<-function(){list(beta=rep(0,n+m),tau.y=1,yf1=rep(0,n+m))}
# inits<-function(){list(alpha=rep(0,n+m),beta=rep(0,n+m),tau.y=1,tau.b=1,tau.a=1,yf1=rep(0,n+m))}
# inits<-function(){list(beta=rep(0,n+m),tau.y=1,tau.b=1,yf1=rep(0,n+m),g=0)}
# inits<-function(){list(beta=rep(0,n),tau.y=1,tau.b=1,yf1=rep(0,n),g=1)}
# inits<-function(){list(beta=rep(0,n),tau.y=1,tau.b=1,yf1=rep(0,n))}
# inits<-function(){list(beta=rep(0,n),yf1=rep(0,n))}

#-Selecting parameters to monitor-
parameters<-c("beta","tau","yf1")
# parameters<-c("beta","tau.y","tau.b","yf1","g")
# parameters<-c("beta","tau.y","tau.b","yf1")
# parameters<-c("alpha","beta","tau.y","tau.b","tau.a","yf1")
# #parameters<-c("beta","tau.y","tau.b","yf1")
# parameters<-c("beta","yf1")

#-Running code-
#OpenBUGS
ej7a.bugs<-bugs(data,inits,parameters,model.file="Ej7a.txt",
               n.iter=50000,n.chains=2,n.burnin=5000)
# ej7a.bugs<-bugs(data,inits,parameters,model.file="Ej7a.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000)
# ej7b.bugs<-bugs(data,inits,parameters,model.file="Ej7b.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000)
# ej7c.bugs<-bugs(data,inits,parameters,model.file="Ej7c.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000,debug=TRUE)
# ej7d.bugs<-bugs(data,inits,parameters,model.file="Ej7d.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000)
#JAGS
ej7a.jags<-jags(data,inits,parameters,model.file="Ej7a.txt",
               n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)

## Warning in jags.model(model.file, data = data, inits = init.values,
## n.chains = n.chains, : Unused variable "m" in data

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 13
##    Unobserved stochastic nodes: 19
##    Total graph size: 173
## 
## Initializing model

# ej7a.jags<-jags(data,inits,parameters,model.file="Ej7a.txt",
#               n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)
# ej7b.jags<-jags(data,inits,parameters,model.file="Ej7b.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)
# ej7c.jags<-jags(data,inits,parameters,model.file="Ej7c.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)

exploracionMCMC(ej7a.bugs)

#Resumen (estimadores)
#OpenBUGS
out.sum<-ej7a.bugs$summary

#JAGS
# out.sum<-ej7a.sim$BUGSoutput$summary

#Tabla resumen
out.sum.t<-out.sum[grep("beta",rownames(out.sum)),c(1,3,7)]
# out.sum.t<-cbind(out.sum.t,apply(out$beta,2,prob))
# dimnames(out.sum.t)[[2]][4]<-"prob"
print(out.sum.t)

##               mean        2.5%     97.5%
## beta[1]  0.2591838 -0.24870000 0.7760050
## beta[2]  1.6631132  0.64189750 2.7300250
## beta[3]  0.2049331 -0.05134025 0.4605025
## beta[4]  2.4290951  1.42100000 3.4850000
## beta[5] -0.4859289 -1.19700000 0.2071000

#DIC
out.dic<-ej7a.bugs$DIC
# out.dic<-ej7a.sim$BUGSoutput$DIC
print(out.dic)

## [1] 3.633

#Predictions
out.yf<-out.sum[grep("yf1",rownames(out.sum)),]
ymin<-min(data$y,out.yf[,c(1,3,7)])
ymax<-max(data$y,out.yf[,c(1,3,7)])
xmin<-min(data$t)
xmax<-max(data$t+m)

#x vs. y
par(mfrow=c(1,1))
plot(data$x,data$y,type="p",col="grey50",ylim=c(ymin,ymax))
points(data$x,out.yf[,1],col=2,pch=16,cex=0.5)
segments(data$x,out.yf[,3],data$x,out.yf[,7],col=2)

#t vs y
par(mfrow=c(1,1))
plot(data$t,data$y,type="b",col="grey80",ylim=c(ymin,ymax),xlim=c(xmin,xmax))
lines(data$t,out.yf[1:n,1],col=2)
lines(data$t,out.yf[1:n,3],col=2,lty=2)
lines(data$t,out.yf[1:n,7],col=2,lty=2)

# lines(data$t[n]:(data$t[n]+m),out.yf[n:(n+m),1],col=4)
# lines(data$t[n]:(data$t[n]+m),out.yf[n:(n+m),3],col=4,lty=2)
# lines(data$t[n]:(data$t[n]+m),out.yf[n:(n+m),7],col=4,lty=2)
#betas
out.beta<-out.sum[grep("beta",rownames(out.sum)),]
ymin<-min(out.beta[,c(1,3,7)])
ymax<-max(out.beta[,c(1,3,7)])
plot(out.beta[,1],type="l",ylim=c(ymin,ymax))
lines(out.beta[,3],lty=2)
lines(out.beta[,7],lty=2)

Modelo B - Dinámico

#-Defining data-
m<-2
# data<-list("n"=n,"m"=m,"y"=milk$y,"x"=milk$x,"t"=milk$t)
# data<-list("n"=n,"m"=m,"y"=milk$y,"x"=milk$x/max(milk$x),"t"=milk$t/max(milk$t))
data<-list("n"=n,"m"=m,"y"=scale(milk$y)[1:n],"x"=scale(milk$x)[1:n],"t"=scale(milk$t)[1:n])
# data<-list("n"=n,"m"=m,"y"=c(scale(milk$y)[1:(n-2)],NA,NA),"x"=scale(milk$x)[1:n],"t"=scale(milk$t)[1:n])

#-Defining inits-
# inits<-function(){list(alpha=0,beta=0,tau=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,5),tau=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,n+m),tau.y=1,tau.b=1,yf1=rep(0,n+m))}
# inits<-function(){list(beta=rep(0,n+m),tau.y=1,yf1=rep(0,n+m))}
# inits<-function(){list(alpha=rep(0,n+m),beta=rep(0,n+m),tau.y=1,tau.b=1,tau.a=1,yf1=rep(0,n+m))}
inits<-function(){list(beta=rep(0,n+m),tau.y=1,tau.b=1,yf1=rep(0,n+m),g=0)}
# inits<-function(){list(beta=rep(0,n),tau.y=1,tau.b=1,yf1=rep(0,n),g=1)}
# inits<-function(){list(beta=rep(0,n),tau.y=1,tau.b=1,yf1=rep(0,n))}
# inits<-function(){list(beta=rep(0,n),yf1=rep(0,n))}

#-Selecting parameters to monitor-
# parameters<-c("beta","tau","yf1")
parameters<-c("beta","tau.y","tau.b","yf1","g")
# parameters<-c("beta","tau.y","tau.b","yf1")
# parameters<-c("alpha","beta","tau.y","tau.b","tau.a","yf1")
# #parameters<-c("beta","tau.y","tau.b","yf1")
# parameters<-c("beta","yf1")

#-Running code-
#OpenBUGS
ej7b.bugs<-bugs(data,inits,parameters,model.file="Ej7b.txt",
               n.iter=50000,n.chains=2,n.burnin=5000)
# ej7c.bugs<-bugs(data,inits,parameters,model.file="Ej7c.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000,debug=TRUE)
# ej7d.bugs<-bugs(data,inits,parameters,model.file="Ej7d.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000)
#JAGS
ej7b.jags<-jags(data,inits,parameters,model.file="Ej7b.txt",
              n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)

## Warning in jags.model(model.file, data = data, inits = init.values,
## n.chains = n.chains, : Unused variable "x" in data

## Warning in jags.model(model.file, data = data, inits = init.values,
## n.chains = n.chains, : Unused variable "t" in data

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 13
##    Unobserved stochastic nodes: 33
##    Total graph size: 70
## 
## Initializing model

# ej7c.jags<-jags(data,inits,parameters,model.file="Ej7c.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)

exploracionMCMC(ej7b.bugs)

dicMCMC(ej7b.bugs)

## [1] -3.792

library(loo)

## This is loo version 2.0.0.
## **NOTE: As of version 2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session. Visit mc-stan.org/loo/news for details on other changes.

## 
## Attaching package: 'loo'

## The following object is masked _by_ '.GlobalEnv':
## 
##     milk

waic(ej7b.bugs$sims.list$yf1)

## Warning: 2 (13.3%) p_waic estimates greater than 0.4. We recommend trying
## loo instead.

## 
## Computed from 90000 by 15 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic      7.3  6.8
## p_waic         3.3  2.0
## waic         -14.6 13.5

## Warning: 2 (13.3%) p_waic estimates greater than 0.4. We recommend trying
## loo instead.

#Resumen (estimadores)
#OpenBUGS
out.sum<-ej7b.bugs$summary

#JAGS
# out.sum<-ej7b.sim$BUGSoutput$summary

#Tabla resumen
out.sum.t<-out.sum[grep("beta",rownames(out.sum)),c(1,3,7)]
# out.sum.t<-cbind(out.sum.t,apply(out$beta,2,prob))
# dimnames(out.sum.t)[[2]][4]<-"prob"
print(out.sum.t)

##                 mean       2.5%      97.5%
## beta[1]  -0.74923274 -1.0410000 -0.3742975
## beta[2]  -0.52620168 -0.8365025 -0.2397000
## beta[3]  -0.35656219 -0.7443000 -0.1099000
## beta[4]  -0.96425303 -1.2160000 -0.5735950
## beta[5]  -0.97862136 -1.2470000 -0.6167975
## beta[6]  -0.96865277 -1.2140000 -0.5414000
## beta[7]  -0.28629377 -0.6011000 -0.0148190
## beta[8]   0.09480327 -0.2546025  0.3489000
## beta[9]  -0.02168113 -0.2781025  0.3265025
## beta[10]  0.28211721  0.0213080  0.6339000
## beta[11]  1.02523909  0.6784000  1.2880000
## beta[12]  1.32703186  1.0430000  1.6510000
## beta[13]  2.15987641  1.6560000  2.4310000
## beta[14]  2.35416292  0.7859900  3.9290000
## beta[15]  2.67084564  0.2567975  5.7800000

#DIC
out.dic<-ej7b.bugs$DIC
# out.dic<-ej7b.sim$BUGSoutput$DIC
print(out.dic)

## [1] -3.792

#Predictions
out.yf<-out.sum[grep("yf1",rownames(out.sum)),]
ymin<-min(data$y,out.yf[,c(1,3,7)])
ymax<-max(data$y,out.yf[,c(1,3,7)])
xmin<-min(data$t)
xmax<-max(data$t+m)

#x vs. y
par(mfrow=c(1,1))
plot(data$x,data$y,type="p",col="grey50",ylim=c(ymin,ymax))
points(data$x,out.yf[1:n,1],col=2,pch=16,cex=0.5)
segments(data$x,out.yf[1:n,3],data$x,out.yf[,7],col=2)

#t vs y
par(mfrow=c(1,1))
plot(data$t,data$y,type="b",col="grey80",ylim=c(ymin,ymax),xlim=c(xmin,xmax))
lines(data$t,out.yf[1:n,1],col=2)
lines(data$t,out.yf[1:n,3],col=2,lty=2)
lines(data$t,out.yf[1:n,7],col=2,lty=2)
lines(data$t[n]:(data$t[n]+m),out.yf[n:(n+m),1],col=4)
lines(data$t[n]:(data$t[n]+m),out.yf[n:(n+m),3],col=4,lty=2)
lines(data$t[n]:(data$t[n]+m),out.yf[n:(n+m),7],col=4,lty=2)

#betas
out.beta<-out.sum[grep("beta",rownames(out.sum)),]
ymin<-min(out.beta[,c(1,3,7)])
ymax<-max(out.beta[,c(1,3,7)])
plot(out.beta[,1],type="l",ylim=c(ymin,ymax))
lines(out.beta[,3],lty=2)
lines(out.beta[,7],lty=2)

Modelo C - Libro Compton

Sin escalar

Utiliza datos sin escalar

#-Defining data-
m<-2
data<-list("n"=n,"m"=m,"y"=milk$y,"x"=milk$x,"t"=milk$t)
# data<-list("n"=n,"m"=m,"y"=milk$y,"x"=milk$x/max(milk$x),"t"=milk$t/max(milk$t))
# data<-list("n"=n,"m"=m,"y"=scale(milk$y)[1:n],"x"=scale(milk$x)[1:n],"t"=scale(milk$t)[1:n])
# data<-list("n"=n,"m"=m,"y"=c(scale(milk$y)[1:(n-2)],NA,NA),"x"=scale(milk$x)[1:n],"t"=scale(milk$t)[1:n])

#-Defining inits-
# inits<-function(){list(alpha=0,beta=0,tau=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,5),tau=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,n+m),tau.y=1,tau.b=1,yf1=rep(0,n+m))}
# inits<-function(){list(beta=rep(0,n+m),tau.y=1,yf1=rep(0,n+m))}
# inits<-function(){list(alpha=rep(0,n+m),beta=rep(0,n+m),tau.y=1,tau.b=1,tau.a=1,yf1=rep(0,n+m))}
# inits<-function(){list(beta=rep(0,n+m),tau.y=1,tau.b=1,yf1=rep(0,n+m),g=0)}
inits<-function(){list(beta=rep(0,n),tau.y=1,tau.b=1,yf1=rep(0,n),g=1)}
# inits<-function(){list(beta=rep(0,n),tau.y=1,tau.b=1,yf1=rep(0,n))}
# inits<-function(){list(beta=rep(0,n),yf1=rep(0,n))}

#-Selecting parameters to monitor-
# parameters<-c("beta","tau","yf1")
parameters<-c("beta","tau.y","tau.b","yf1","g")
# parameters<-c("beta","tau.y","tau.b","yf1")
# parameters<-c("alpha","beta","tau.y","tau.b","tau.a","yf1")
# #parameters<-c("beta","tau.y","tau.b","yf1")
# parameters<-c("beta","yf1")

#-Running code-
#OpenBUGS
# ej7a.bugs<-bugs(data,inits,parameters,model.file="Ej7a.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000)
# ej7a.bugs<-bugs(data,inits,parameters,model.file="Ej7a.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000)
# ej7b.bugs<-bugs(data,inits,parameters,model.file="Ej7b.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000)
ej7c.bugs<-bugs(data,inits,parameters,model.file="Ej7c.txt",
               n.iter=50000,n.chains=2,n.burnin=5000)
# ej7d.bugs<-bugs(data,inits,parameters,model.file="Ej7d.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000)
#JAGS
# ej7a.jags<-jags(data,inits,parameters,model.file="Ej7a.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)
# ej7a.jags<-jags(data,inits,parameters,model.file="Ej7a.txt",
#               n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)
# ej7b.jags<-jags(data,inits,parameters,model.file="Ej7b.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)
ej7c.jags<-jags(data,inits,parameters,model.file="Ej7c.txt",
               n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)

## Warning in jags.model(model.file, data = data, inits = init.values,
## n.chains = n.chains, : Unused variable "m" in data

## Warning in jags.model(model.file, data = data, inits = init.values,
## n.chains = n.chains, : Unused variable "t" in data

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 13
##    Unobserved stochastic nodes: 29
##    Total graph size: 89
## 
## Initializing model

exploracionMCMC(ej7c.bugs)

dicMCMC(ej7c.bugs)

## [1] 41.18

#Resumen (estimadores)
#OpenBUGS
out.sum<-ej7c.bugs$summary

#JAGS
# out.sum<-ej7a.sim$BUGSoutput$summary

#Tabla resumen
out.sum.t<-out.sum[grep("beta",rownames(out.sum)),c(1,3,7)]
# out.sum.t<-cbind(out.sum.t,apply(out$beta,2,prob))
# dimnames(out.sum.t)[[2]][4]<-"prob"
print(out.sum.t)

##               mean   2.5%  97.5%
## beta[1]   9.753333  9.619  9.892
## beta[2]  10.040285  9.897 10.160
## beta[3]  10.231874 10.070 10.350
## beta[4]  10.152901 10.030 10.310
## beta[5]  10.319832 10.190 10.460
## beta[6]  10.430720 10.310 10.600
## beta[7]  10.904846 10.750 11.030
## beta[8]  11.142070 10.990 11.270
## beta[9]  11.264554 11.140 11.420
## beta[10] 11.537864 11.400 11.680
## beta[11] 11.873193 11.710 12.000
## beta[12] 11.953647 11.830 12.120
## beta[13] 12.332323 12.170 12.470

#DIC
out.dic<-ej7c.bugs$DIC
# out.dic<-ej7a.sim$BUGSoutput$DIC
print(out.dic)

## [1] 41.18

#Predictions
out.yf<-out.sum[grep("yf1",rownames(out.sum)),]
ymin<-min(data$y,out.yf[,c(1,3,7)])
ymax<-max(data$y,out.yf[,c(1,3,7)])
xmin<-min(data$t)
xmax<-max(data$t+m)

#x vs. y
par(mfrow=c(1,1))
plot(data$x,data$y,type="p",col="grey50",ylim=c(ymin,ymax))
points(data$x,out.yf[,1],col=2,pch=16,cex=0.5)
segments(data$x,out.yf[,3],data$x,out.yf[,7],col=2)

#t vs y
par(mfrow=c(1,1))
plot(data$t,data$y,type="b",col="grey80",ylim=c(ymin,ymax),xlim=c(xmin,xmax))
lines(data$t,out.yf[1:n,1],col=2)
lines(data$t,out.yf[1:n,3],col=2,lty=2)
lines(data$t,out.yf[1:n,7],col=2,lty=2)
lines(data$t[n]:(data$t[n]+m),out.yf[n:(n+m)],col=4)
lines(data$t[n]:(data$t[n]+m),out.yf[n:(n+m)],col=4,lty=2)
lines(data$t[n]:(data$t[n]+m),out.yf[n:(n+m)],col=4,lty=2)

#betas
out.beta<-out.sum[grep("beta",rownames(out.sum)),]
ymin<-min(out.beta[,c(1,3,7)])
ymax<-max(out.beta[,c(1,3,7)])
plot(out.beta[,1],type="l",ylim=c(ymin,ymax))
lines(out.beta[,3],lty=2)
lines(out.beta[,7],lty=2)

Escalados

Pésima predicción

#-Defining data-
m<-2
# data<-list("n"=n,"m"=m,"y"=milk$y,"x"=milk$x,"t"=milk$t)
# data<-list("n"=n,"m"=m,"y"=milk$y,"x"=milk$x/max(milk$x),"t"=milk$t/max(milk$t))
data<-list("n"=n,"m"=m,"y"=scale(milk$y)[1:n],"x"=scale(milk$x)[1:n],"t"=scale(milk$t)[1:n])
# data<-list("n"=n,"m"=m,"y"=c(scale(milk$y)[1:(n-2)],NA,NA),"x"=scale(milk$x)[1:n],"t"=scale(milk$t)[1:n])

#-Defining inits-
# inits<-function(){list(alpha=0,beta=0,tau=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,5),tau=1,yf1=rep(1,n))}
# inits<-function(){list(beta=rep(0,n+m),tau.y=1,tau.b=1,yf1=rep(0,n+m))}
# inits<-function(){list(beta=rep(0,n+m),tau.y=1,yf1=rep(0,n+m))}
# inits<-function(){list(alpha=rep(0,n+m),beta=rep(0,n+m),tau.y=1,tau.b=1,tau.a=1,yf1=rep(0,n+m))}
# inits<-function(){list(beta=rep(0,n+m),tau.y=1,tau.b=1,yf1=rep(0,n+m),g=0)}
inits<-function(){list(beta=rep(0,n),tau.y=1,tau.b=1,yf1=rep(0,n),g=1)}
# inits<-function(){list(beta=rep(0,n),tau.y=1,tau.b=1,yf1=rep(0,n))}
# inits<-function(){list(beta=rep(0,n),yf1=rep(0,n))}

#-Selecting parameters to monitor-
# parameters<-c("beta","tau","yf1")
parameters<-c("beta","tau.y","tau.b","yf1","g")
# parameters<-c("beta","tau.y","tau.b","yf1")
# parameters<-c("alpha","beta","tau.y","tau.b","tau.a","yf1")
# #parameters<-c("beta","tau.y","tau.b","yf1")
# parameters<-c("beta","yf1")

#-Running code-
#OpenBUGS
# ej7a.bugs<-bugs(data,inits,parameters,model.file="Ej7a.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000)
# ej7a.bugs<-bugs(data,inits,parameters,model.file="Ej7a.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000)
# ej7b.bugs<-bugs(data,inits,parameters,model.file="Ej7b.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000)
ej7c.bugs<-bugs(data,inits,parameters,model.file="Ej7c.txt",
               n.iter=50000,n.chains=2,n.burnin=5000)
# ej7d.bugs<-bugs(data,inits,parameters,model.file="Ej7d.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000)
#JAGS
# ej7a.jags<-jags(data,inits,parameters,model.file="Ej7a.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)
# ej7a.jags<-jags(data,inits,parameters,model.file="Ej7a.txt",
#               n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)
# ej7b.jags<-jags(data,inits,parameters,model.file="Ej7b.txt",
#                n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)
ej7c.jags<-jags(data,inits,parameters,model.file="Ej7c.txt",
               n.iter=50000,n.chains=2,n.burnin=5000,n.thin=1)

## Warning in jags.model(model.file, data = data, inits = init.values,
## n.chains = n.chains, : Unused variable "m" in data

## Warning in jags.model(model.file, data = data, inits = init.values,
## n.chains = n.chains, : Unused variable "t" in data

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 13
##    Unobserved stochastic nodes: 29
##    Total graph size: 89
## 
## Initializing model

exploracionMCMC(ej7c.bugs)

dicMCMC(ej7c.bugs)

## [1] 38.04

#Resumen (estimadores)
#OpenBUGS
out.sum<-ej7c.bugs$summary

#JAGS
# out.sum<-ej7a.sim$BUGSoutput$summary

#Tabla resumen
out.sum.t<-out.sum[grep("beta",rownames(out.sum)),c(1,3,7)]
# out.sum.t<-cbind(out.sum.t,apply(out$beta,2,prob))
# dimnames(out.sum.t)[[2]][4]<-"prob"
print(out.sum.t)

##                 mean   2.5%     97.5%
## beta[1]  -0.35127338 -1.225 0.4876025
## beta[2]  -0.25270435 -1.073 0.5766000
## beta[3]  -0.22759774 -1.140 0.6793025
## beta[4]  -0.35547787 -1.817 0.8075025
## beta[5]  -0.17950499 -1.711 1.3060000
## beta[6]  -0.01181980 -1.420 1.6980250
## beta[7]  -0.04091683 -1.339 1.4390250
## beta[8]  -0.13862157 -1.400 1.1830000
## beta[9]  -0.16708131 -1.238 0.9510050
## beta[10] -0.34236650 -1.367 0.7015025
## beta[11] -0.74512861 -2.004 0.4988075
## beta[12] -1.08827754 -3.001 0.6342025
## beta[13] -1.62006552 -5.055 0.5313025

#DIC
out.dic<-ej7c.bugs$DIC
# out.dic<-ej7a.sim$BUGSoutput$DIC
print(out.dic)

## [1] 38.04

#Predictions
out.yf<-out.sum[grep("yf1",rownames(out.sum)),]
ymin<-min(data$y,out.yf[,c(1,3,7)])
ymax<-max(data$y,out.yf[,c(1,3,7)])
xmin<-min(data$t)
xmax<-max(data$t+m)

#x vs. y
par(mfrow=c(1,1))
plot(data$x,data$y,type="p",col="grey50",ylim=c(ymin,ymax))
points(data$x,out.yf[,1],col=2,pch=16,cex=0.5)
segments(data$x,out.yf[,3],data$x,out.yf[,7],col=2)

#t vs y
par(mfrow=c(1,1))
plot(data$t,data$y,type="b",col="grey80",ylim=c(ymin,ymax),xlim=c(xmin,xmax))
lines(data$t,out.yf[1:n,1],col=2)
lines(data$t,out.yf[1:n,3],col=2,lty=2)
lines(data$t,out.yf[1:n,7],col=2,lty=2)
lines(data$t[n]:(data$t[n]+m),out.yf[n:(n+m)],col=4)
lines(data$t[n]:(data$t[n]+m),out.yf[n:(n+m)],col=4,lty=2)
lines(data$t[n]:(data$t[n]+m),out.yf[n:(n+m)],col=4,lty=2)

#betas
out.beta<-out.sum[grep("beta",rownames(out.sum)),]
ymin<-min(out.beta[,c(1,3,7)])
ymax<-max(out.beta[,c(1,3,7)])
plot(out.beta[,1],type="l",ylim=c(ymin,ymax))
lines(out.beta[,3],lty=2)
lines(out.beta[,7],lty=2)

Sesión 15

Ejemplo 7

EDA

Modelo O

Modelo A

Modelo B - Dinámico

Modelo C - Libro Compton

Sin escalar

Escalados